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ABSTRACT 


Accurate Modeling of flexible space structures is an important field that is 
currently under investigation. Parameter estimation, using methods such as max- 
imum likelihood, is one of the ways that the model can be improved. The maximum 
likelihood estimator has been used to extract stability and control derivatives 
from flight data for many years. Most of the literature on aircraft estimation 
concentrates on new developments and applications, assuming familiarity with 
basic estimation concepts. This paper presents some of these basic concepts. 

The paper briefly discusses the maximum likelihood estimator and the aircraft 
equations of motion that the estimator uses. The basic concepts of minimization 
and estimation are examined for a simple computed aircraft example. The cost 
functions that are to be minimized during estimation are defined and discussed. 
Graphic representations of the cost functions are given to help illustrate the 
minimization process. Finally, the basic concepts are generalized, and estima- 
tion from flight data is discussed. Specific examples of estimation of struc- 
tural dynamics are included. Some of the major conclusions for the computed 
example are also developed for the analysis of flight data. 


INTRODUCTION 


Accurate modeling of flexible space structures is an important area that is 
currently under investigation. The mathematical nr -deling of these structures 
can be improved using parameter estimation. Such techniques have been success- 
fully used to estimate aircraft stability and control derivatives and refine 
aircraft mathematical models. Some of the experience gained in the aircraft 
problem can be applied directly to analysis of flexible space structures. 


The maximum likelihood estimator has been used to obtain stability and con- 
trol estimates from flight data for nearly 20 years. The results of many appli- 
cations have been reported worldwide. Reference 1 contains a representative 
list of some of these reports. Several good texts (including Refs. 2 and 3) 
contain thorough treatments of the theory of maximum likelihood estimation. 
Experience reports (Refs. 1, 4, and 5) pointing out practical considerations for 
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applying the maximum likelihood estimator have also been published* Stability 
and control derivatives estimated f roa flight data are currently required for 
correlation studies with predictive techniques, handling qualities documentation, 
design compliance, aircraft simulator enhancement and refinement, and control 
system design. Correlation, simulation, and control system design applications 
(including the space shuttle) are discussed in Ref. 6. Current studies have 
concentrated on estimation model structure determination (Refs. 7 and 8), equa- 
tion error with state reconstruction (Refs. 9 to 11), and maximum likelihood 
estimation in the frequency domain (Refs. 12 and 13). 

Most of the reports in the estimation area concentrate on new developments 
and applications, assuming familiarity with the basic concepts of maximum like- 
lihood estimation. In this paper some of these basic concepts are reviewed, 
concentrating on simple, idealized models. These simple models provide insights 
applicable to a wide variety of real problems. 

This paper also presents sooe of the basics of maximum likelihood estimation 
as applied to the aircraft problem. It briefly discusses the maximum likelihood 
estimator and the aircraft equations of notion that the sstinator uses. The 
basic aspects of minimization and estimation are then examined in detail for a 
simple computed aircraft example. Finally, the discussion is expanded to the 
general aircraft estimation problem including specific examples of estimation of 
structural dynamics. 


SYMBOLS 


A,B|C,D, F ,G 


a N 


a x 

a y 

a z 

b 

c* 

C n 

Cx 

C Y 

C Z 


system matrices 

normal acceleration positive upward, g 
longitudinal acceleration, g 
lateral acceleration, g 
normal acceleration positive upward, g 
reference span, ft 
coefficient of rolling moment 
coefficient of yawing moment 
coefficient of axial force 
coefficient of side force 
coefficient of normal force 


f(*), g(-) 


general functions 



GG* 

g 

H 

lx* *xz 

i 

a 

Kp 

L 


L* 


l YJ 

M 


N 


n 

a 

Pg 

p 

q 

q 

R 

r 

s 

T 

t 

u 

V 


measurement noise covariance matrix 

acceleration due to gravity, ft/sec 2 

approximation to the information matrix 

moment of inertia about subscripted axis, slug-ft 2 

general index 

cost function 

sidewash factor 

rolling moment divided by X x , deg/sec 2 
rolling moment, ft-lb 
rolling moment due to yaw jet, ft-lb 
pitching moment divided by Iy, deg/sec 2 
mass, slug 

number of time points or cases or yawing moment divided 
by I z , deg/sec 2 

state noise vector or number of unknowns 

estimated roll rate due to turbulence, deg/sec 

roll rate, deg/sec 

pitch rate, deg/sec 

dynamic pressure, lb/ft 2 

innovation covariance matrix 

yaw : ce, g/oec 

reference area, ft 2 

time increment, sec 

time, sec 

control input vector 
forward velocity, ft/sec 
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X 


x ay*yay» 2 ay 


H 

a 

e 

A 

A 

& 

«a 

«e 

fir 

n 

e 

u 

c 

a 

T 

♦ 

;■ 

0) 

Subscripts : 

p,q,r,a,a,0,fi, 

fi*fia*fir*fi« 


0 

m 


state vector 

distance between lateral accelerometer and the 
center ol gravity along the appropriate axis, ft 

observation vector 

predicted Kalman-filtered estimate 

angle of attack, deg 

angle of sideslip, deg 

estimated angle of sideslip due to turbulence, deg 

time sample interval, sec . 

control deflection, deg 

aileron deflection, deg 

elevon deflection, deg 

rudder deflection, deg 

measurement noise vector 

pitch angle, deg 

mean 

vector of unknowns 
standard deviation 
time, sec 

transition matrix or bank angle, deg 

integral of transition matrix, or heading angle, deg 

Irequency, rad/sec 

partial derivative with respect to subscripted quantity 

bias or at time zero 
measured quantity 
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Other nomenclature: 


predicted estimate 

estimate 

transpose 

indicates moment in ft- lb 


MAXIMUM LIKELIHOOD ESTIMATION 


The concept of maximum likelihood is discussed in this section. First the 
general heuristic problem is discussed, and then the specific equations for 
obtaining maximum likelihood estimates for the aircraft problem are given. In 
the following sections, both the concepts and the computations involved in a 
simple but realistic example are discussed in detail. 

The aircraft parameter estimation problem can be defined quite simply in 
general terms. The system investigated is assumed to be modeled by a set of 
dynamic equations containing unknown parameters. To determine the values of the 
unknown parameters, the system is excited by a suitable input, and the input and 
actual system response are measured. The values of the unknown parameters are 
then inferred based on the requirement that the model response to the given 
input match the actual system response. When formulated in this manner, the 
problem of identifying the unknown parameters can be easily solved by many 
methods; however, complicating factors arise when application to a real system 
is considered. 

The first complication results from the impossibility of obtaining perfect 
measurements of the response of any real system. The inevitable sensor errors 
are usually included as additive measurement noise in the dynamic model. Once 
this noise is introduced, the theoretical nature of the problem changes drasti- 
cally. It is no longer possible tc exactly identify the values of the unknown 
parameters; instead, the values must be estimated by some statistical criterion. 
The theory of estimation in the presence of measurement noise is relatively 
straightforward for a system with discrete time observations, requiring only 
basic probability. 

The second complication of real systems is the presence of state noise. 

State noise is random excitation of the system from unmeasured sources, the 
standard example for the aircraft stability and control problem being 
atmospheric turbulence. If state noise is present and measurement noise is 
neglected, the analysis results in the regression algorithm. 

When both state and measurement noise are considered, the problem is more 
complex than in the cases that have only state noise or only measurement noise. 
Reference 14 develops a mixed continuous/discrete maximum likelihood formulation 
that allows for both state and measurement noise. This formulation has a con- 
tinuous system model with discrete sampled observations. 
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The final problem for real systems is modeling* It has been assumed through- 
out the above discussion that for some value (called the "correct" value) of the 
unknown parameter vector, the system is correctly described by the dynaadc model* 
Physical systems are seldom described exactly by single dynamic models, so the 
question of modeling error arises. No comprehensive theory of modeling error is 
available. The most common approach is to ignore it: Any modeling error is 

simply treated as state noise or measurement noise, or both, in spite of the 
fact that the modeling error may be deterministic rather than random. The 
assumed noise statistics can then be adjusted to include the contribution of 
the modeling error. This procedure is not rigorously justifiable, but, combined 
with a carefully chosen model, it is probably the best approach available. 

With the above discussion in mind, it is possible to make a more precise, 
mathematically probabilistic statement of the parameter estimation problem. The 
first step is to define the general system model (aircraft equations of motion). 
This model can be written in the continuous/discrete form as 

x(to) * xo (1) 

i(t) - f[x(t),u(t),£] + F(£)n(t) (2) 

z(ti) - g[x(ti),u(ti),ei + GU)hi (3) 

where x is the state vector, z is the observation vector, f and g are system 
state and observation functions, u is the known control input vector, £ is the 
unknown parameter vector, n is the state noise vector, and f) is the measurement 
noise vector. The state noise vector is assumed to be zero-mean white Gaussian 
and stationary, and the measurement noise vector is assumed to be a sequence of 
independent Gaussian randou variables with zero mean and identity covariance. 

For each possible estimate of the unknown parameters, a probability that the 
aircraft response time histories attain values near the observed values can then 
be defined. The maximum likelihood estimates are defined as those that maximize 
this probability. Maximum likelihood estimation has many desirable statistical 
characteristics: for example, it yields asymptotically unbiased, consistent, and 
efficient estimates (Ref. 15). 

If there is no state noise and the matrix G is known, then the maximum 
likelihood estimator minimizes the cost function 

1 N ~ , 

J(£) - t 2 [sE(ti) - z^(t^) ]*(GG*)”1 (z(ti) - z^tti)] (4) 

i«1 

where GG* is the measurement noise covariance matrix, and z^(ti) is the computed 
response estimate of z at t^ for a given value of the unknown parameter vector 

£. The cos t function is a function of the difference between the measured and 
computed time histories. 
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If Eqs. (2) and (3) are linearizsd (as is the case for the stability and 
control derivatives in the aircraft problem), 


x(to) 

- XQ 

(5) 

x(t) 

* Ax(t) + Bu( t) + Fn( t) 

(€) 

z(ti) 

■ Cx(ti) + Dutt^) + GTii 

(7) 


For the* no-state-noise case, 

x$(t 0 ) 
*e< ti+i> 

h (t i> 


the term of Bq. (4) can be approximated by 


x 0 («) 


(8) 

♦«5(ti) 

+ 'l'tu(ti) + u(t i+1 )]/2 

(9) 

Cx^tt*) 

+ Du(ti) 

(10) 


where 


♦ 

* 


exp [A(t i+1 - t ± ) ] 



exp (At) dT o 


When state noise is important, the nonlinear form of Eqs. (1) to (3) is 
intractable. For the linear model defined by Eqs. (5) to (7), the cost function 
that accounts for state noise is 


J(U * 1 E lz(ti) - Z£(t< ))*R- 1 lz(ti) - zr(ti)] + 1 N In |r| (11) 

i«1 2 

where R is the innovation covariance matrix# The Z£(t^) term in Bq# (11) 
is the Kalman-filtered estimate of z, which f if the state noise covariance 
is zero, reduces to the form of Bq. (4)# If there is no state noise, 
the second term of Eq. (11) is of no consequence (unless one wishes to include 
elements of the G matrix as unknowns), and R can be replaced by GG* which makes 
Eq. (11) the same as Eq. (4). 

To minimize the cost function J(£)# we can apply the Hewton-Raphson 
algorithm which chooses successive estimates of the vector of unknown coef- 

A A 

ficients, £. Let L be the iteration number. The L + 1 estimate of 5 is then 
obtained from the L estimate as follows: 

£l+i - £l - tv|j(£ L ) )-i [v|j(i L ) ] d2) 
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The first and second gradients are defined as 


N ^ 

V £J(S> = - 2 [z(ti) - zsttiHMGG*)- 1 [V^z € ( t A ) ] (13) 

i*1 

2 N 

^J(C) * 2 [V 5 z C (t i )]*(GG*)- 1 [V 5 2 ? (t i )l 
i=1 

N 

- 2 [z(ti) - z e (ti)]MGG*)-1[vf zgtti)] (14a) 

i=1 


The Gauss-Newton approximation to the second gradient is 

2 N 

V<£) 2 J IV^z^ttiJJ^tGG^J-’lV^z^tti)] 
i=1 


(14b) 


The Gauss-Newton approximation, which is sometimes referred to as modified 
Newton-Raphson, is computationally much easier than the Newton-Raphson approxi- 
mation because the second gradient of the innovation never needs to be calcu- 
lated* In addition, it can have the advantage of speeding the convergence of 
the algorithm, as is discussed in the SIMPLE AIRCRAFT EXAMPLE section. 

Figure 1 illustrates the maximum likelihood estimation concept. The meas- 
ured response of the aircraft is compared with the estimated response, and the 
•.ifference between these responses is called the response error. The cost func- 
tions of Eqs. (4) and (11) include this response error. The Gauss-Newton com- 
putational algorithm is used to find the coefficient values that maximize the 
cost function. Each iteration of this algorithm provides a new estimate of the 
unknown coefficients on the basis of the response error. These new estimates of 
the coefficients are then used to update the mathematical model of the aircraft, 
providing a new estimated response and, therefore, a new response error. The 
updating of the mathematical model continues iterati/ely until a convergence 
criterion is satisfied. The estimates resulting from this procedure are the 
maximum likelihood estimates. 

The maximum likelihood estimator also provides a measure of the reliability 
of each estimate based on the information obtained from each dynamic maneuver. 
This measure of the reliability, analogous to the standard deviation, is called 
the Cramer-Rao bound (Ref. 16) or the uncertainty level. The Cramer-Rao bound 
as computed by current programs should generally be used as a measure of rela- 
tive accuracy rather than absolute accuracy. The bound is obtained from the 
approximation of the information matrix, H. This matrix equals the approxima- 
tion to the second gradient given by Eg. (14b). The bound for each unknown is 
the square root of the corresponding diagonal element of H. That is, for 

the ith unknown, the Cramer-Rao bound is /H(i,i) . 
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The Maine-Iliff formulation (Ref, 14) and minimization algorithm discussed 
above are implemented with the Iliff-Maine code (MMLE3 maximum likelihood esti- 
mation program). The program and computational algorithms are described fully 
in Ref, 17, All the computations shown and described in the remainder of the 
paper use the algo-rithms exactly as described in Ref. 17. 


AIRCRAFT EQUATIONS OF MOTION 


For the discussion tnat follows in later sections of this paper, some ’ owl- 
edge of the airjraft equations of motion is assumed. To clarify some of tnc*i 
discussion, the aircraft equations are discussed briefly in this section. 

First, the axis system on which the aircraft equations of motion are based 
xs discuss id* Figure 2(a) shows thj aircraft reference xxly-axis system and the 
conventional control r^rfaces. The origin of the body-axis system is at the 
center of gravity. The sign convention for this axis system is defined by the 
right-hand rule with the x-axis defined as positive forward on the aircraft* 

The longitudinal acceleration (a x ) and nondimensional axial force coefficient 

(Cy) are defined along this axis, and the roll rate (p) and rolling moment (L 1 ) 

are defined about this axis* The y-axis is defined as positive out the right 

wing. The lateral acceleration (ay) and nondime'isional side force coefficient 

(Cy) are defined along this axis, and the pitch rate (q) and pitching moment 

(M' ) are defined about this axis. The z-axis is defined as positive out the 

bottom of the aircraft. The normal acceleration (a 2 ) and nondimensional normal 
force coefficient (Cg) are defined along this axis, and the yaw rate (r) and 

yawing moment (N* ) are defined about this axis* The normal acceleration is 
sometimes defined is positive upward but is then referred to as a^« The three 
moments (L*, M 1 , and N') are usually divided by the corresponding moments of 
inertia (I x , Iy, and I 2 }, and are then referred to without the prime as L, M, 

and N. These quantities are nondimensionalized (C 4 , C m , and C n , respectively) 

for use in tb- equations of motion soon to be discussed. The primary control 
about t;v- roll axis (x-axis) is the aileron ( 6 a ), about the pitch axis (y-axis) 

is the ele/ator ( 6 e ), and about the yaw axis (z-axis) is the rudder ( 6 r ). Some 
aircraft have other controls, but in this paper these will only be defined where 
they are discussed (the reaction control jets on the space shuttle, for example). 

The Euler angles ♦, 0, and + define the aircraft attitude with respect to 
the earth. These angles define the rotations which transform earth-fixed axes 
to the aircraft reference body-axis system of Fig* 2(a)* . The order of rotation 
must be about the z-axis (♦), then the y-axis ( 6 ), and finally the x-axis ($) 
for the aircraft equations of motion that will be written subsequently* 

For stability and control analysis, the velocity of the aircraft with 
respect to the air (not with respect to the earth) is of primary interest* 

Figure 2(b) shows the relationship between the aircraft axis system and the flow 
angles* The flow angle in the x-z plane is the angle of attack (a), and the 
flow angle in the x-y plane is the angle of sideslip (3). A more rigorous and 
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detailed definition is required for the derivation of the equations of motion, 
but the above definitions are sufficient to define the following equation of 
motion. 

Generalized nonlinear equations of motion are given in detail in Ref. 17, 
which fully describes the Iliff-Haine code (MMLE3 program). All computations 
an' 1 aircraft examples in this paper use the linearized form for the lateral- 
directional equations. These equations are given below and referred to in the 
remainder of the paper. 


where 


where the 6 



+ 0 O ) 


COS 

0 sin ♦ + p 

sin 

a - r cos a 

(15) 

pl x - 

• 

• rl^ 

* qsbCji 

+ qr(Iy - I 

z) + 

M 1 xz 

(16) 

riz - 

e 

■ P*xz 

* qsbc n 

♦ pqu x - * 

y - 

q rI xz 

(17) 

• 

♦ - p 

+ r cos ♦ 

tan 

0 + q sin t 

tan 

8 + ^0 

(18) 

2y - 

c v g 8 - 

>c *i 


+ C *r TV + 

C Y{ 5 

+ c Yq 

(19) 

c i * Cjtg3 

+ C *P 

fl 

+ C i: 

r ^ + C *5* 

6b 

+ c* 0 + 

(20) 

c n “ c ng£ 

+ Cn P 


+ C n, 

r "3V + Cn S 6 

+ c r. 0 ♦ c ng 

(21) 


term is summed over all controls. 


The observation equations are 


zg X0 

* Kg 0 - Tj— P + y— r 
Pm * P 


» r 


♦m * ♦ 

. y .aic, -4 s 

Ym mg g 


. *a y . 

p * -T' 


^ (.2 . r*) 


• • • 

Pm * P + PO 

• • * 

* r + *0 


(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 
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The state, control, and observation vectors for the lateral-directional node 
can then be defined as 


X = (8 p r <fr)* 

(29) 

u = (6 a 6 r )* 

(30) 

• • 

z * < e m Pm J'm ♦* *y m Pm r m>* 

(31) 


SIMPLE AIRCRAFT EXAMPLE 


The basic concepts involved in a parameter estimation problem can be illus- 
trated by using a simple example representative of a realistic aircraft problem. 
The example chosen here is representative of an aircraft that exhibits pure 
rolling motion from an aileron input. Hiis example, although simplified, typi- 
fies the motion exhibited by many aircraft in particular flight regimes, such as 
the F-14 aircraft flying at high dynamic pressure, the F-111 aircraft at moder- 
ate speeds with the wing in the forward position, and the T-37 aircraft at low 
speed. 

Derivation of an equation describing this motion is straightforward. 

Figure 2(c) shows a sketch of an aircraft with the x-axis perpendicular to the 
plane cf the figure (positive forward on the aircraft). The rolling moment (L'), 
roll rate (p), and aileron deflection (6 a ) are positive as shown. For this 

example, the only state is p and the only control is 6 a » The result of summing 
moments is 

IxP-L'(p,fia) (32) 

The first-order Taylor expansion then becomes 

p - Lpp + L$ a « a (33) 

where 

L'« I X L 

Since the aileron is the only control, it is notationally simpler to use 6 
Instead of 6 a for the discussion of this example. Equation (33) can then be 
written as 

P - Lpp + l$6 (34) 

; alternate approach that results in the same equation is to combine Eq. (16) 
with £q. (20), substituting for C;, end then eliminate the terms that are zero 
for our example. This yields 

pl x = gab Cf p + Cj^e (35) 
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where p is the roll rate and 6 is the aileron deflection. Rearranging terms# 
the equation can be put into the dimensional derivative form of En. (34). 

Equation (34) is a simple aircraft equation where the forcing function is 
provided by the aileron and the damping by the damping-in-roll term, Ip. In 

subsequent sections we examine in detail the parameter estimation problem where 
Bq. (34) describes the system. For this single-degree -of -freedom problem, the 
maximum likelihood estimator is used to estimate either Lp or ox both for a 
given computed time history* 

We will assume that the system has measurement noise, but no state noise as 
in Eqs. (1), (2), and (3). Equation (4) then gives the cost function for maximum 
likelihood estimation. The weigh cing QG* is unisqportant for this problem, so 
let it equal 1. For our example, Eqs. (2) and (3) become x^ * p^ and z± * x^. 
Therefore, Bq. (4) becomes 


J(Ip,L<$) 



fpi - pi(V l 5 )j 2 


(36) 


where pi is the value of the measured response p at time t^ and Pj.(Ip,I>6) is the 

a a 

computed time history of p at time t^ for Ip * Ip and Throughout the 

rest of the paper, where computed data (not flight data) are used, the measured 

time history refers to p^, and the computed time history refers to Pi(Ip,L()« 

The computed time history is a function of the current estimates of lp and L$, 
but the measured tine history is not. 


The most straightforward method of obtaining pi is with Eqs. (3) and (8). 

In terms of the notation stated above, 

Pi+1 * ♦Pi ♦ ♦ (*! + *i+1>/ 2 (37) 


where 


♦ = exp (IpA) 


f A 1*6(1 - e*P (WO! 

♦ * I exp (LpT) dx • ; * 

J o r 

and A is the length of the sample interval (t^+i - t^). Simplifying the 
notation 


then 


*1+1/2 * r *i + *1+1 )/2 
Pi+1 = ♦Pi + ^ 1 + 1/2 


(38) 


(39) 
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The aaxinum likelihood estimate is obtained by ainimizing Eq. (36). The 
Gauss-Newton method described earlier is used for this miniaization. Equa- 
tion (12) is used to determine successive values of the estimates of the 
unknowns during the minimization. 

a A A A A 

Por this simple problem, £ * [Lp L$]* and successive estimates of Lp e.id 

are determined by updating Bq. (12). The first and second gradients of Bq. (12) 
are defined by Bqs. (13) and (14). The complete set of equations is given in 
Ref. 17. 

The entire procedure can now be written for obtaining the maximum likelihood 
estimates for this simple example. To start the algorithm, an initial estimate 

A 

of Lp and L$ is needed. This is the value of Co* With Bq. (12), Ci end sub- 

sequently Cl are defined by using the first and second gradients of J(Lp,L£) 

from Eq. (36). The gradients for this particular example from Eq. (13) and 
'14b) are 

„ N ~ 

V£J(£l> * -2 ( Pi - Pi>*£ Pi (40) 

i*1 

N 

V^J(i L ) H 2 (41) 

i*1 

With the specific equations defined in this section for this simple example, 
we can now proceed in the next section to the computational details of a speci- 
fic example. 


Computational Details of Minimization 

In the previous section we specified the equations for a simple example and 
described the procedure for obtaining estimates of the unknowns from a dynamic 
maneuver* In this section we give the computational details for obtaining the 
estimates. Some of the basic concepts of parameter estimation are best shown 
with computed data where the correct answers are known. Therefore, in this sec- 
tion we study two examples involving computed time histories. The first example 
is based on data that have no measurement noise, which results in estimates that 
are the same as the correct value. The second example contains significant 
measurement noise; consequently, the estimates are not the same as the correct 
values. Throughout the rest of the paper, wLere computed data are used, the 
term "no-noise case" is used for the case with no noise added and "noisy case" 
for the case where noise has been added. 

Since we are studying a simple computed example, it is desirable to keep it 
simple enough to complete some or all of the calculations on a home computer or, 
wit* some labor, on a calculator. With this in mind, the number of data points 
netds to be kept sn^ll. Por this computed example, 10 points (time samples) are 
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used. Tl.e simulated data, which we refer to as the Measured data, are based on 
Eq. (34). Me use the same correct values of Lp and Lj (-0.2500 and 10.0, 
respectively) for both examples. In addition, the same input (6) is used for 
both e> Moles . the saaple interval (A) is 0.2 sec, and the initial condi tiona 
are zeic- Tables of all the significant interaediate values are given with each 
example. These values are given to four significant digits, although to obtain 
exactly tl saae values with a computer or calculator requires the use of 13 
signified- digits, as in the computation of these tables. If the four-digit 
nua- er t are used in the computation, the answers will be a few tenths of a 
percent off, but will still serve to illustrate the minimization accuracy. In 

A 

both examples, the initial values of Lp and L$ (or are -0e5 and 15.0, 
respectively# 

Example With No Measurement Noise 

The leasarement time history for no measurement noise (no-noise case) is 
shown in Fig. 3# The aileron input starts at zero, goes to a fixed value, and 
then returns to zero. The resulting roxl-cate time history is also shown. Hie 
values of the measured roll rate to 13 sigrificant digits are given in Table 1 
along with the aileron input. 

A A 

Table 2 shows the values for Lp, L$, and J for each iteration, along with 
the values of A a..d 4> needed for calculations of p£. in three iterations the 
algorithm converges to the correct values to four significant digits for both Lp 

A 

and L$. L$ overshoots slightly on the first iteration and then cones quickly to 

at 

the correct answer. Lp overshoots slightly on the second iteration. 

Figure 4 shows the match between the measured data and the computed data for 
each of the first three terations. The match is very good after two iterations. 
The match is nearly exact after three iterations. 

Although the .a lgorithm has converged to four-digit accuracy in Lp and L$, 

the value of the cost function, J, continues to decrease rapidly between itera- 
tions 3 and 4. This is a consequence of using the maximum likelihood estimator 
on data with no measurement noise. Theoretically, using infinite accuracy the 
value of J at the minimum should le zero. However, with finite accuracy the 
value of J becomes small but never quite zero. This value is a function of the 
number of significant digits that are being used. For the 13-digit accuracy 

o o 

used here, the cost eventually decreases to approximately 0.3 * 10~*. 

Example With Measurement Noise 

The data used _n this example (noisy case) are the same as those used in the 
previous section, except that pseudo-Gaussian noise has been added to the roll 
rate. The time history i*. a .iorn in Fig. 5. The signal-to-noise ratio is quite 
low in thi~ example, as is readily apparent by comparing Figs. 3 and 5. The 

exact values of the time history to 13-digit accuracy are shown in ^able 3. The 

A A 

values of Lp, L£, $, 1>, and J are shown for each iteration in Table 4. The 
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algorithm converges in four iterations. The behavior of the coefficients as 
they approach convergence is much like the no-noise case. The most notable 

A A 

results of this case are the converged values of Lp and L$, which are somewhat 
different from the correct values. The match between the measured and computed 
time h: tc-y i* shown in Fig. 6 for each iteration. No change in the match is 
apparent for the last two iterations. The match is very good considering the 
amount of measurement noise. 

In Fig. 7, the computed time history lor the correct values of Lp and Lg is 
compared to that for the noisy-case estimates of T and Lg. Because the 
algorithm converged to values somewhat different w an the correct values, the 
two compute . time histories are similar but not identical. 

The accuracy of the converged elements can be assessed by looking at the 
Cramer-Rao inequality (Refs. 16 and 17) discussed earlier. The Cramer-Rao bound 
can be obtained from the following approximation to the information matrix. 

H * 2(J, iniaum )(v|j>-V(H-1) 


The Cramer-Rao bounds for Lp and Lg are the square roots of the diagonal ele- 
ments of the H matrix, or /H(1,1) and /h( 2,2), respectively. The Cramer-Rao 

A A 

bounds are 0*1593 and 1.116 for Ip and L$, respectively. The errors in Lp and 
Lfi are less than the bounds. 


Cost Functions 


In the previous section we obtained the maximum likelihood estimates for 
computed time histories by minimizing the values of the cost function. To fully 
understand what occurs in this minimization, we must study in more detail the 
form of the cost functions a.id some of their more important characteristics. In 
this section, the cost function for the no-noise case is discussed briefly. The 
cost function of the noisy case is then discussed in more detail. The same two 
time histories studied in the previous section are s-amined here. The noisy 
case is more interesting because it has a meaningful Cramer-Rao bound and is 
more representative of aircraft flight data. 

First we will look at the one-dimensional case where L$ is fixed at the 

correct value, because it is easier to grasp some of the characteristics of the 
cost function in one dimension. Than we will look at the two-dimensional case, 
where both Lp and Lj are varying. It is important to remember that everything 

shown in this paper on cost functions is based on computed time histories that 
are defined by Eq. (36). For every time history we might choose (computed or 
flight data), a complete cost function is defined. For the case of n variables, 
the cost function defines a hypersurface of n + 1 dimensions. It might occur to 
us that we could just construct this surface and look for the minimum, avoiding 
the need to bother with the minimization algorithm* This is not a reasonable 
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approach because, in general, the number of variables is greater than two* 
Therefore, the cost function can be described mathematically but not pictured 
graphically. 


One-Dimensional Case 

To illustrate the many interesting aspects of cost functions, it is eas- 
iest to first look at cost functions having one variable. In an earlier sec- 
tion, the cost function of Lp and Lg was minimized* That cost function is most 

interesting in the Lp direction. Ther^Jore, the one-variable cost function 
studied here is J(Lp). All subsequent discussions are for J(Lp) with Lg equal 

to the correct value of 10. Figure 8 shows the cost function plotted as a func- 
tion of Lp for the case where there is no measurement noise (no-noise case). As 

expected for this case, the minimum cost is zero and occurs at the correct value 
of Lp s -0.2500. It is apparent that the cost increases much more slowly for a 

more negative Lp than for a positive Lp. In fact, the slope of the curve tends 
to become less negative where Lp is more negative than -1.0. Physically this 
makes sense since the more negative values of Lp represent cases of high damping, 
and the positive Lp represents an unstable system. Therefore, the Pi for posi- 
tive Lp becomes increasingly different from the measured time history for small 
positive increments in Lp. For very large damping (very negative Lp) the system 

would show essentially no response. Therefore, large increases in damping 
result in relatively small changes in the value of J(Lp). 

In Fig. 9, the cost function based on the time history with measurement 
noise (noisy case) is plotted as a function of Lp. The correct value of 
Lp (-0.2500) and the value of Lp (-0.3218) at the minimum of the cost (3.335) 

are both indicated on the figure. The general shape of the cost function in 
Fig. 9 is similar to that shown in Fig. 8. Figure 10 shows the comparison 
between the cost functions based on the time histories with and without measure- 
ment noise. The comments relating to the cost function of the no-noist. ^ase 
also apply to the cost function based on the noisy case. Figure 10 shows 
clearly that the two cost functions are shifted by the difference the value 
of Lp at the minimum and increased by the difference in the minimum cost. One 

would expect only a small difference in the value of the cost when far from the 
minimum. This is because the "estimated" time history is so far from the 
measured time history that it becomes irrelevent as to whether the measured time 
history has noise added* Therefore, for large values of cost, the difference in 
the two cost functions should be small in comparison to the total cost. 

Figure 11 shows the gradient of J(Lp) plotted as a function of Lp for the 

noisy case. This is the function for which we were trying to find the zero (or 
equivalently, the minimum of the cost function) using the Gauss -Newton method of 
a previous section. The gradient is zero at L p * -0.3218, which corresponds to 
the value of the minimum of J(Lp). 

The difference between the Newton-Raphson method (Eq. (14a)) and the 
Gauss-Newton method (Eq. (14b)) of minimization has been mentioned previously. 
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For this simple one -dimensional case, we can easily compute the second gradient 
both with the second term of Eq. (14a) (Newton-Raphson), and without the second 
term (Gauss -Newton, Eq. (14b)). Figure 12 shows a comparison between the 
Newton-Raphson and the Gauss -Newton approximation second gradients. The 
Gauss-Newton second gradient (dashed line) always remains positive because it 
is the sum of quadratic terms (squared for the one-dimensional example). The 
Newton-Raphson second gradient can be positive or negative, depending upon the 
value of the second partial with respect to Lp. Other than the difference in 
sign for the more negative Lp, the two curves have similar shapes. 

As stated earlier, the Gauss-Newton method can be shown to be superior to 
Newton-Raphson in certain cases. We can demonstrate obvious cases of this 
with our example. An easy way to select a spot where problems with the Newton- 
Raphson method will occur is to look for places where the second gradient (slope 
of the gradient) is near zero o^ negative. Figure 11 has such a region near 
Lp = -1.0. If we choose a point where the gradient slope is exactly zero, we 

are forced to divide by zero in Eq. (12) with the Newton-Raphson meti ,jd. This 
point is at Lp = -1.13 in Fig. 12. If the value of the slope of the gradient is 

negative, then the Newton-Raphson method will go to very negative values of Lp. 
For very negative values of Lp, the cost becomes asymptotically constant and the 
gradient becomes nearly zero. In that region, the Newton-Raphson algorithm 
would diverge towards negative infinity. If the slope of the gradient is 
positive but small, we still have a problem with the Newton-Raphson me < hod. 

Figure 13 shows the first iteration starting from Lp « -0.95 for both Gauss- 

Newton and Newton-Raphson. The Newton-Raphson method selects a point where the 
tangent of the gradient at Lp * -0.95 intersects the zero line. This results in 

the selection of an Lp of approximately 2.6 in the first iteration. From that 
value it requires many iterations to return to the actual minimum. On the other 
hand, the Gauss-Newton method selects a value for Lp of approximately -0.09 and 

converges to the minimum to four-digit accuracy in two more iterations. With 
more complex examples a comparison of the convergence properties of the two 
algorithms becomes more difficult to visualize, but the problems are generaliza- 
tions of the situation we have observed with the one-dimensional example. 

The usefulness of the Cram&r-Rao bound was discussed in the Example With 
Measurement Noise section. At this point it is useful to digress briefly to 
discuss some of the ramifications of the Cramer-Rao bound for the one-dimensional 
case. The Cramer-Rao bound only has meaning for the noisy case. In the noisy 
example, the estimate of Lp is -0.3218 and the Cramer-Rao bound is 0.0579. The 

calculation of the Cramer-Rao bound was defined in the previous section for both 
one-dimensional and two-dimensional examples. The Cramer-Rao bound is an esti- 
mate of the standard deviation of the estimate. One would expect the scatter in 
the estimates of Lp to be of about the same magnitude as the estimate of the 

standard deviation. For the one-dimensional case discussed here, the range 
(Lp (-0.3218) plus or minus the Cramer-Rao bound (0.0579)) nearly includes the 

correct value of Lp (-0.2500). If noisy cases are generated for many time his- 
tories (adding different measurement noise to each time history), then the sam- 
ple mean and sample standard deviation of the estimates for these cases can be 
calculated. Table 5 gives the sample mean, sample standard deviation, and the 
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standard deviation of the sample mean (standard deviation divided by the square 
root of the number of cases) for 5, 10, and 20 cases* The sample mean, as 

expected, gets closer to the correct value of -0*2500 as the number of cases 
increases. This is also reflected by the decreasing values in column 4 of 
Table 5, which are estimates of the error in the sample mean* Column 3 of 
Table 5 shows the sample standard deviations, which indicate the approximate 
accuracy of the individual estimates* This standard deviation, which stays mora 
or less constant, is approximately equal to tne Cramer-Rao bound for the noisy 
case being studied here. In fact, the Cramer-Rao bounds for each of the 20 
noisy cases used here (not shown in the table) do not change much from the 
values found for the noisy case being studied. Both of these results are in 
good agreement with the theoretical characteristics (Ref. 16) of the Cramer-Rao 
bounds and maximum likelih od estimators in general. 

The examples shown here indicate the value of obtaining more sample time 
histories (maneuvers). More samples improve confidence in the estimate of the 
unknowns. The same result holds true in analyzing actual flight time histories 
(maneuvers); thus it is always advisable to obtain several maneuvers at a given 
flight condition to improve the best estimate of each derivative. 

The size of the Cramer-Rao bounds and of the error between the correct value 
and the estimated value of Lp is determined to a large extent by the length of 
the time history and the amount of noise added to the correct time history. For 
the example being studied here, it is apparent from Fig. 5 that the amount of 
noise being added to the time history is large. The effect of the power of the 

A 

measurement noise (GG*, Eqs. (3) and (4)) on the estimate of Lp (that is, Lp) 
for the time history is given in Table 6. The estimate of Lp is much improved 

by decreasing the measurement noise power. A reduction in the value of G to 
one-tenth of the valu^ in the noisy example being studied yields an acceptable 
estimate of Lp. For l light data, the measurement noise is reduced by improving 
the accuracy of the output of the measurement sensors. 

Two-Dimensional Case 

In this section the cost function (which is deiendent on both l p and i»6 ) is 
studied. The no-noise case is examined first, followed by the noisy case. 

No-noise case . Even though the cost function is a function of only two 
unknowns, it is much more difficult to visualize than the one-unknown case. The 
cost function over a reasonable range of Lp and L$ is shown in Fig. 14. The 

cost increases very rapidly in the region of positive Lp and large values of 
Lj. The reason is just an extension of the argument fcr positive Lp given in 
the previous section. The shape of the surface can be depicted in greater 
detail if we examine only the values of the cost function less than 200 for Lp 

less than 1.0. Figure 15 shows a view of this restricted surface from the upper 
end of the surface. The minimum must lie in the curving valley that gets 
broader as we go to the far side of the surface. Now that we have a picture of 
the surface, we can look at the isoclines of constant cost on the Lp-versus-Lj 

plane. These isoclines are shown in Fig. 16. The minimum of the cost function 
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is inside the closed isocline. The steepness of the cost function in the 
poeitive-Lp direction is once again apparent. Inside the closed isocline the 

shape is more nearly elliptical, indicating chat the cost is nearly quadratic 
here, so fairly rapid convergence in this region would be expected. The Lp axis 

becomes an asymptote in cost as L 5 approaches zero. The cost is constant for 
L5 - 0 because no response would result from any aileron input. The estimated 

response is zero for all values of Lp, resulting in constant cost. 

Figure 16 shows the region of the minimum value of the cost function, which, 
as seen in the earlier example (Table 2), occurs at the correct values for Lp 
and L$ of -0.2500 and 10, respectively. This is also evident by looking at the 

cost function surface shown in Fig. 17. The surface has its minimum at the 
correct value. As expected, the value of the cost function at the minimum is 
zero. 

Noisy case . As shown before in the one-dimensional case, the primary dif- 
ference between the cost functions for the no-noise and noisy cases was a shift 
in the cost function. In that instance, the noisy case was shifted so that the 
minimum was at a higher cost and a more negative value of Lp. In the two- 

dimensional case, the no-noise and noisy cost functions exhibit a similar shift. 
For two dimensions the shift is in both the Lp and L$ directions. The shift is 

small enough that the difference between the two cost functions is not visible 
at the scale shown in Fig. 14 or from the perspective of Fig. 15. Figure 18 
shows the isoclines of constant cost for the noisy case. The figure looks much 
like the isoclines for the no-noise case shown in Fig. 16. The difference 
between Figs. 16 and 18 is a shift in Lp of about 0.1. This is the difference 

in the value of Lp at the minimum for the no-noise and noisy cases. Heuristi- 

cally, one can see that the same would be true for cases with more than two 
unknowns. The primary difference between the two cost functions is near the 
minimum. 

The next logical part of the cost function to examine is near the minimum. 
Figure 19 shows the same view of the cost function for the noisy case as was 
shown in Fig. 17 for the no-noise case. The shape is roughly the same as that 

shown in Fig. 17, but the surface is shifted such that its minimum lies over 

Lp ■= -0.3540 and L$ * 10.24, and is shifted upward to a cost function value of 

approximately 3.3. 

To get a more precise idea of the cost of the noisy case near the minimum, 

we once again need to examine the isoclines. The isoclines (Fig. 20) in this 

region are much more like ellipses than they are ir Fig*>- 16 and 18. We can 
follow the path of the minimization example used before by including the results 
from Table 4 on Fi^. 20. The first iteration (L * 1 ) brought the values of Lp 
and 1*5 very close to the values at the minimum. The next iteration essentially 

selected the values at the minimum when viewed at this scale. One of the rea- 
sons the convergence is so rapid in this region is that the isoclines are nearly 
elliptical, demonstrating that the cost is very nearly quadratic in this region. 
If we had started the Gauss-Newton algorithm at a point where the isoclines 
are tfuch less elliptical (as in some of the border regions in Fig. 18), the 
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convergence would have been much slower initially, but much the same as it 
entered the nearly quadratic region of the cost function. 


Before concluding our examination of the two-dimensional case, we need to 
examine the Cramer-Rao bound. Figure 21 shows the uncertainty ellipsoid, which 

is based on the Cramer-Rao bounds defined in an earlier section. The relation- 
ships between the Cramer-Rao bound and the uncertainty ellipsoid are discuss ;d 
in Ref. 16. The uncertainty ellipsoid almost includes the correct value of Lp 

and Ii< 5 . The Cramer-Rao bound for Lp and L$ can be determined from the projec- 
tion of the uncertainty ellipsoid onto the Lp and L 5 axes, and compared with the 
values given earlier, which were 0.1593 and 1.116 for Lp and L$, respectively. 

ESTIMATION USING FLIGHT DATA 


In the previous several sections we examined the basic mechanics of obtain- 
ing maximum likelihood estimates from computed examples with one or two unknown 
parameters. Now that we have a grasp of these basics, we can explore the esti- 
mation of stability and control derivatives from actual flight data. For the 
computationally much more difficult situation usually encountered using actual- 
flight data, we will obtain the maximum likelihood estimates with th« Ilirf- 
Maine code (MMLE3 program) described in Ref. 17, The equations of motion that 
are of interest are given in the AIRCRAFT EQUATIONS OF MOTION section of this 
paper; the remainder of the equations are given in Ref. 17. 

In general, flight data estimation is fairly complex, and codes as the 

Iliff-Maine code must usually be used to assist in the analysis* However, one 
must still be cautious about accepting the results; that is# the estimates must 
fit the phenomenology, and the match between the measured and computed time his- 
tories must be acceptable. This is true in all flight regimes, but one must be 
particularly careful in potential problem situations such as ( 1 ) in separated 
flow at high Mach numbers or high angle of attack, (2) with unusual aircraft 
configurations such as the oblique wing (Ref. 18), or (3) with modern high- 
performance aircraft with high-gain feedback loops. In any of the above cases, 
one should be particularly careful where there are even small anomalies in the 
match. These anomalies may indicate ignored terms in the equations of motion, 
separated flow, nonlinearities, sensor problems, insufficient resolution 
(Ref. 1), sensor location (Ref. 1), time or phase lags (Refs. 1 and *9), or aay 
of a long list of ocher problems. 

The following brief examples are intended to show how the above caveats and 
the computed examples of previous sections can be used to assist in the analy- 
sis. Ii. the computed example, the desirability of low-noise sensors, an ade- 
quate model, and several maneuvers at a given flight condition is shown. 


Hand Calculation Example 

Sometimes evaluation of a fairly complex flight maneuver can be augmented 
with a simple hand calculation. One example of this can be found for the space 
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shuttle. The space shuttle is a large double-delta-winged vehicle designed to 
enter the atmosphere from space and land horizontally. The entry control system 
consists of 12 vertical reaction-control-systen (RCS) jets (six up-firing and 
six down-firing), 8 horizontal RCS jets (four left-firing and four right-firing) 
4 elevon surfaces, a body flap, and a split rudder surface. The locations of 
these devices are shown in Pig. 22. the vertical jets and the elevens are used 
for both pitch and roll control. The jets and elevons are used symae tries lly 
for pitch control and asymmetrically for roll control. The space shuttle con- 
trol system is described briefly in Ref. 6. 

The shuttle example used here is from a maneuver obtained at a Mach number 
of approximately 21 and an angle of attach of approximately 40°. The controls 
being used for this lateral-directional maneuver are the differential elevons 
and the side-firing jets (yaw jets). The maneuver is shown in Fig. 23. 

Equations (15) to (31) describe the equations of motion. A simplified approach 
can be used to determine some of the derivatives by hand. The approach is one 
that has been used since the beginning of dynamic analysis of flight maneuvers. 
In particular, for this maneuver the slope of the rates can be used to determine 
the yaw jet control derivatives. This is possible for this example, even with a 
high-gain feedback system, because the yaw jets are essentially step functions, 
and the slope of the rates p and r can be determined before the vehicle and the 
differential elevon (aileron) responses become significant. The rolling moment 
due to yaw jet (Lyj) is particularly important for the shuttle (Ref. 6 discusses 
the essential nature of flight-determined Lyj in the redefinition of entry 

maneuvers) and is, in general, more difficult to obtain than the more dominant 
yawing moment due to yaw jet. Therefore, as an illustrative example, Lyj is 

determined by hand. Figure 24 shows yaw jet activity and smoothed roll rate 
plotted at expanded scales. The equation for Lyj is given by 

Lyj ■ pl x / (Number of yaw jets) (42) 

$ - Ap/At - ♦ (0.1 ) (43) 


Therefore, given that I x 2 900,000 slug-ft 2 , and the number of yaw jets is 4, 

L yj 2 2750 ft-lb. 

The same maneuver was analyzed with MMLE3, and the resulting match 4 shown 
in Fig. 25# The match is very good except for a small mismatch in p at about 
6 sec. This small mismatch was studied separately with MMLE3 and found to be 
caused by a nonlinearity in the aileron derivative. The value from MMLE3 for 
Lyj is 2690 ft-lb, which for the accuracy used here is essentially the same 

value as obtained by the simplified method. The aileron derivatives would be 
difficult to determine as accurately as the yaw jet derivatives. Although good 
estimates can seldom be obtained with the slope method discussed here, rough 
estimates can usually be obtained to gain some insight into values obtained with 
MMLE3 (or any other maximum likelihood program). These rough estimates can then 
be used to help explain unexpected values of estimates from an estimation 
program. 
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Sometimes a flight example becomes too complex to allow anything other than 
qualitative estimates to be determined by hand. The example shown in Fig. 26 is 
the determination of the rudder derivative for tl e F-8 aircraft with the yaw 
augmentation system on. This example, taken from Ref. 20, includes an aileron 

pulse and a rudder pulse. Although an independent pilot rudder pulse is input 
during the maneuver, the rudder is largely responding to the lateral accelera- 
tion feedback. When the rudder is moving, several other variables are also 
moving, thus making it difficult to use the simplified approach just discussed. 
However, C n ^ can be roughly determined when the rudder moves, approximately 

1.7 sec from the start of the maneuver. Most of the slope of yaw rate 4 8 pro- 
bably caused by the rudder, but a poor estimate would be obtained using the hand 
calculation. 


Cost Function for Full Aircraft Problem 


The analysis of a lateral-directional maneuver obtained in flight typically 
has from 15 to 25 unknown parameters (as shown in Eqs. (15) and (31)), in 
contrast to the one or two in the simple aircraft example. This makes detailed 
examples unwieldy and any graphic presentation of the cost function impossible. 
Therefore, in this section we are primarily examining tl,e estimation procedure 
and the process of the minimization. 

For our flight example, we have chosen a lateral-directional maneuver, with 
both aileron and rudder inputs, that has 17 unknown parameters. The data ar<« 
from the oblique wing aircraft (Ref. 18) with the wing unskewed during the 
maneuver. This example was chosen because it is a typical maneuver. The time 
history of the data and the subsequent output of MMLE3 have been published in 
Ref. 21. Some results of the analysis are shown in Table 7. The match between 
the measured time history (solid ^ines) and the estimated (calculated) time 
history (dashed lines) is shown as a function of iteration in Fig. 27. Fig- 
ures 27(a) to (e) are for iterations 0 to 4, respectively. Table 7 shows that 
the cost remains unchanged after four iterations. A similar result was obtained 
for the two-dimensional simple aircraft example in Fig. 6 and Table 4. 

Of the many things the analyst must consider in obtaining estimates, the two 
most important or.sg are how good is the match and how good is the convergence. 

A satisfactory match and monotor convergence are necessary, but not suffi- 
cient, conditions for a success.' x analysie- Figure 27(e), although not per- 
fect, is a very good match. The convergence can best be evaluated by looking at 
the normalized cost in the last row of Table 7. The cost has converged rapidly 
and monotonically in four iterations, and it remains at the converged cost. 

These factors are convincing evidence that the convergence is complete. There- 
fore. the criteria of match and convergence are satisfied in our example. In 
some cases we might encounter cost that does not converge rapidly (in four to 
six iterations) or monotonically, or stay "exactly" at the minimum value. These 
situations usually indicate at least a small problem in the analysis. These 
problems, if found, are usually traced to an instrumentation or data aquisition 
problem, an inadequate mathematical model, or a maneuver that contains a mar- 
ginal amount of information. 
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Table 7 also shows that the startup values of all the coefficients are zero 
for the control and bias variables* Wind tunnel estimates could have been used 
for starting value ^ but the convergence of the algorithm is not very dependent 

or the startup values. As part of the startup algorithm, the MMLE3 program nor- 
mally holds the derivatives of the state variables constant until after the 
first iteration, is evident in Table 7, 

Figure 27(a) shows the match between the measured and computed data for the 
startup values. The match is very poor because the startup values for the con- 
trol derivatives are all zero, so the only motion is in response to the initial 
conditions. The control derivatives and blase , are determined on the first 
iteration, resulting in the much improved match shown in Fig. 27(b). The match 
after two iterations, shown in Fig. 27(c), is improved as the program further 
modifies the control derivatives and, for the first time, adjusts the deriva- 
tives affecting the natural frequency (C n g and C$g). By the third iteration 

(Fig. 27(d)), the improvement in the match is. almost complete, because minor 
adjustments to the frequency are made and the damping derivatives are changed. 
Fig. 27(e) shows the match when all but the most minor derivatives have ceased 
to change. 

Several general observations can be made based on this well behaved example. 
The strong or most important coefficients hare essentially converged in thr^e 
iterations. The same effect was seen in the simple example — that is, con- 
verged faster than Lp (Table 4). Some of the less important or second-order 
coefficients have only converged to two places after three iterations and are 
still changing by one digit in the fourth place at the end of six iterations. 
Another observation is that for some coefficients (0$^, C n< j , and Cj^ ) even 

though the sign is wrong after the first iteration, the algorithm quickly 
selects their correct values once the important derivatives have stabilized. 

In general, if the analysis of a maneuver has gone well, we do not need to 
spend much time inspecting a table analogous to Table 7. However, i/ there have 
been problems in convergence or in the quality of the fit, a detailed inspection 
of such a table may be necessary. The data may show an important coefficient 
going unstable at an early iteration, which could cause problems later. If ths 
starting values are grossly in error, the algorithm is driven a long way from 
reasonable values and then for many reasons does not behave well. Occasionally 
the algorithm alternately selects from two diverse set of values of two or more 
coefficients on successive iterations, behaving as if the shape of the cost 
function were a narrow multidimensional valley analogous to but more extreme 
than the two-dimensional valley shown in Figs. 18 and 20. 


Cramer-Rao Bounds 


The earlier sections regarding the computed example have shown that the 
Cramer-Rao bound is a good indicator of the accuracy of an estimated parameter. 
The Cramer-Rao bounds can be used in a similar, but somewhat more qualitative, 
fashion on flight data. The Cramer-Rao bounds that are included in MMLE3 (as 
well as many other maximum likelihood estimation programs) have been useful in 
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determining whether estimates are good or bad. The aircraft example discussed 
here has b ^n reported previously (for example, in Refs. 1 and 16). However, 

this example of the use of the Cramer-Rao bound in the assessment of flight- 
derived estimates is pertinent to the thrust of this paper. Figure 28 shows 
estimates of C n ^ as a function of angle of attack for the PA-30 twin-engine 

general aviation aircraft (Ref. 22) at three flap settings. There is a signifi- 
cant amount of scatter, which makes the reliability of the information on C n ^ 

questionable. The data shown are the estimates from the MMLE3 program, which 
also provides the Cramer-Rao bounds for each estimate. Past experience (Ref. 1) 
has shown that if the Cramer-Rao bound is multiplied by a scale factor (the 
result sometimes being called the uncertainty le^el (Refs. 1 and 16)) and plot- 
ted as a vertical bar with the associated estimate, it helps in the interpreta- 
tion of flight-determined results. Figure 29 shows the same data as Fig. 28, 
with the uncertainty levels now included as vertical bars. The estimates with 
small uncertainty levels (Cramer-Rao bounds) are the best estimates, as was 
discussed earlier in the section on Cramer-Rao bounds for the one-dimensional 
case. The fairing shown in Fig. 29 goes through the estimates with small 
Cramer-Rao bounds and ignores the estimates with large bounds. One can have 
great confidence in the fairing of the estimates, because the fairing is well 
defined and consistent when the Cramer-Rao bound information is included. In 
this particular instance, the estimates with small bounds were from maneuvers 
where the aileron forced the motion, and the large bounds were from maneuvers 
where the rudder forced the motion. Therefore, in addition to aiding in the 
fairing of the estimates, the Cramer-Rao bounds help show that the aileron- 
forced maneuvers are superior for estimating C n ^ for the PA-30 aircraft. 

This example illut> crates that the Cramer-Rao bounds are a useful tool in 
assessing flight-determined estimates, just as they were found useful for the 
simple aircraft example with computed datae 


Atmospheric Turbulence ( S tate Noise) 

Atmospheric turbulence (state noise) cannot always be avoided in flight; 
therefore, it is desirable to be aole to obtain stability and control deriva- 
tives in the presence of turbulence. In addition, an estimate of the turbulence 
time history can be of interest, particularly in the implementation of tur- 
bulence suppression systems. 

Many years ago it was demonstrated that the stability and control deriva- 
tives can be adequately determined with maximum likelihood estimation techniques 
for maneuvers performed in smooth air. If the^e techniques, which do not 
account for turbulence, are applied to data obtained in turbulence, net only are 
the resulting matches of the time histories unsatisfac tory but the estimated 
coefficients are unacceptable (Refs. 23 to 25). The technique described in 
Refs. 14, 23, and 25 can account for the effect of turbulence. With this tech- 
nique, maximum likelihood estimates of the stability and control derivatives as 
well as estimates of the turbulence time histories are >t>tained by minimizing 
the cost function given by Fq. (11). Results of the application of the tech- 
nique to longitudinal maneuvers obtained in turbulence have been reported 
previously (Refs. 23 to 25). 



Tne lateral-directional equations (Eqs. (15), (16), (17), (18), and (29)) 
can be modified in a Banner similar to that used to modify the longitudinal 
equations in Refs. 23 to 25. The turbulence (state noise) model is the Dryden 
expression, which is described in Ref. 26. The Iliff-Maine code (Ref. 17) can 
be used to obtain the maximum likelihood estimates where state noise is present. 

Thirty-eight seconds of data from the PA-30 aircraft flying in turbulence 
was analyzed at 50 samples/sec. The best match that could be obtained with the 
maximum likelihood estimation method that does not account fcr turbulence is 
shown in Fig. 30. The match is unacceptable and resulted in poor estimates of 
the stability and control derivatives. Figure 31 shows the match obtained with 
the maximum likelihood estimation technique that accounts for turbulence 
(Re' .14 and 17). The match is excellent and the maneuver provided acceptable 
es mated stability and control derivatives. It is also of interest to compare 
the power spectra of the estimated turbulence time histories. The power 

spectrum of the turbulence component affecting angle of sideslip, shown 

in Fig. 32. Figure 33 presents the power spectrum of the turbulence component 
affecting roll rate, Pg. The slopes of the asymptotes shown in Figs. 32 and 33 
are those defined by the Dtyden expression given in Ref. 26. Good agreement is 

a * 

showit between the power spectra and the asymptotes for 8g and pg. 

The algorithm used here is based on a linearized system described by 
Eqs. (5) to (7) and solved by minimizing the cost function given by Bq. (11). 

The system need not resemble that for the aircraft stability and control problem 
other than in the requirement for linearity. Therefore, many formulations for 
the structural problem are written in the form of Bqs. (5) to (7), and the 
algorithm under discussion can be directly applied with these formulations* 


ESTIMATION FOR SIMPLE STRUCTURAL PROBLEM 


The problem of the flexible space structure is most fully characterized as a 
distributed parameter system with its associated distributed system control 
laws. The model will vary depending upon changes in its configuration or its 
environment, such as solar heating. As in most cases, the preferred solution is 
the simplest successful approach. The lumped system approach is much simpler 
and computationally far more efficient chan the fully distributed parameter sys- 
tem approach. For example, structural mode control base^ on current state-of- 
the-art approaches has proved very successful* Admittedly, the aircraft struc- 
ture is heavier than most spacecraft, but many aircraft structures are highly 
complex, consisting of many subetructures within the main structure* To the 
novice, many of the sp^ce structures currently Seing investigated appear simpler 
than modern, large aircraft. If the lumped parameter system approach used for 
rhe aircraft pro* am is found to be inadequate, it seems likely that distributed 
parameter estimation codes will evolve to whatever complexity is necessary to 
sol Vi the flexible space structure problem. 
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This paper has discussed some of the erqterience gained f roe the application 
of aircraft stability and control analysis to flight data. The codes used for 
this analyses are for lumped paraeeter systems in the time domain. The codes 

have been used successfully for structural problems and are fully adaptable to 
the frequency domain if that is found to be preferable. 

Although few results have been obtained for time-domain structural analysis 
at the Ames Dryden Flight Research Facility, some superficial experience in 
structural time-domain analysis has been obtained. The following two examples 
show how the techniques being used for stability and control analysis can be 
applied to simple structural problems. The preceeding section discussed the 
incorporation of state noise in the model. The following examples do not 
include the use of state noise, but state noise, if warranted, could easily be 
incorporated in the types of examples to be discussed. 


Estimation of Structural Characteristics 


All aircraft have observable structural modes. These modes usually cause no 
difficulty in estimating stability and control derivatives because the struc- 
tural frequencies are higher than the aerodynamic frequencies. In general, if 
the structural frequencies are higher than the highest aerodynamic frequency by 
more than a factor of 5 to 10, they can be neglected unless their amplitude is 
so large as to mask measurement? desired for the aerodynamic analysis. However, 
if one or more structural modes are affecting the aerodynamic modes, as may 
occur in large aircraft, these structural nodes must be included in the mathe- 
matical model being analyzed. 

Even though no completely satisfactory practical results are available that 
account for structural modes and their interactions with the aerodynamics , it is 
interesting to assess the time-domain maximum likelihood analysis of the struc- 
tural modes independent of any interaction. This can be done where a structural 
mode is observed and no significant coupling is apparent. 

Figure 34 shows a structural mode on the lateral acceleration of an aircraft 
where little effect was observed for structural-aerodynamic coupling. The fre- 
quency of the mode is high enough that the mode does not interact with the aero- 
dynamic modes. Therefore, the stability and control derivatives were obtained 
separately and held constant for the succeeding analysis. The analysis con- 
sisted of using the maximum likelihood estimation program KMLE 3 (Ref. 17) with 
a sixth-order model that included the lateral-directional aerodynamic modes plus 
one structural mode. The dynamic pressure and the velocity were allowed to vary 
in the analysis. The structural mode frequency and damping were estimated as 
linear functions of dynamic pressure. The initial conditions were also esti- 
mated. A structural mode frequency of 7.84 Hz was chosen to start the estima- 
tor process. The comparison between the original data and the match obtained 
wiu.\ the maximum likelihood estimation method is shown in Fig. 35. The two time 
histories are in good agreement at the beginning of the maneuver and at the end 
of the maneuver, but they are 180° out of phase at a time of approximately 
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0*3 sec. The match shown in Fig* 35 suggests that the maximum likelihood esti- 
mator has reached a local minimum but not the global minimum. Multiple minima 
are not normally a problem when obtaining the stability and control derivatives 
of aircraft with the maximum likelihood estimation method. 

The reason for the multiple minima is demonstrated by the following simple 
scalar example. Let the noiseless measured response be z(t) * sin (u)Qt) the 
estimated response be Z£ * sin (u)t), where w is the only unknown coefficient. 
Then, by Eq. (4), the cost function becomes 

j(u>,t) = f [sin (u)Qt) - sin (wt)]^ dt 
J 0 

= * " 4^ sin (2 "0 t) " i sin <2u> T ) 


2u) w 0 

- — sin (wT) cos (wnT) - cos (wT) sir (unT) 

0)2 _ wq2 w 

If T is chosen to represent 10 cycles, as shown in Fig* 35, then for an u>q of 
1 rad/sec, T equals 20w. in Fig. 36, the cost function J(<i>,20t) is shown as a 
function of a). The global minimum is at an o> of 1 rad/sec, as it should be, but 
there are many local minima at increments of approximately 0*05 rad/sec. If a 
value of less than 0.97 or greater than 1*03 were chosen for a starting estimate 
of (i), the algorithm would converge to a local minimum* If a value of between 
0*98 and 1.02 were chosen, it would converge to the global minimum* Therefore, 
for this example where 10 cycles were observed, the starting value of u> must be 
less than 3 percent from the correct anwser to converge to the global minimum* 

Figure 37 shows a sine wave for the global minimum along with a sine wave 
with a frequency that varies 10 percent from the global minimum. The sine waves 
are in phase at the beginning and end, ana 180° out of phase in the middle* 

These data appear similar to those shown for flight data in Fig. 35. If only 
one or two cycles were used for the analysis, the problem illustrated in Fig. 37 
would be minimized* This is apparent in Fig. 38 where only the first cycle of 
Fig. 37 is shown. 

If T is chosen to represent only one cycle and wq remains equal to 1 rad/sec 

(as in Fig. 38), then T equals 2x* The cost function J(w,2*) is presented as a 

function of u> in Fig. 39* The global minimum is correctly at an w of 1 rad/sec, 

but now the algorithm converges to the global minimum if u> is started within 

approximately 25 percent of the correct value. 

Knowing the sensitivity of the algorithm when a record with many lightly 
damped cycles is being analyzed, the data of Fig. 34 can be reanalyzed starting 
closer to the observed frequency. Starting the maximum likelihood estimation 
method with an u) of 9.0 results in the fit shown in Fig. 40. This is an accep- 
table fit of the data. 


223 



Based on the preceding results, if data are to be analyzed where many cycles 
of a structural node are present, the structural mode frequency, u>, must be clo- 
sely approximated before starting the estimation process. 


Structural Modes in Space 

In the process ^f analyzing aircraft flight data, the authors have fre- 
quently observed results that clearly exhibit unmodeled dynamics. The unmodeled 
dynamics could be caused by many phenomena, such as higher-order aerodynamic 
modes or structural modes. These modes can usually be ignored and left unmo- 
deled because they have no effect on the results of primary interest in tne ana- 
lysis. If the unmodeled modes cannot be ignored, then the system equations must 
be revised to include the unmodeled modes. 

The authors have not yet found it necessary to model structural modes for 
data obtained in space in the process of obtaining control derivatives for the 
space shuttle. However, the structural modes have been observed. Figure 41 
shows the response of the space shuttle to the firing of a roll jet and a yaw 
jet at an altitude of 430,000 ft. The space shuttle configuration and the loca- 
tion of the RCS jets are shown in Fig. 22. The changes in the rigid-body rates 
and lateral acceleration caused by the jet firings are apparent in Fig. 41. The 
structural modes are also excited by the jets, as evidenced by the increased 
ringing in each signal at the time of the jet firings. The roll jet firing has 
little effect on the rigid-body response for the yaw rate and lateral accelera- 
tion; nowever, the yaw jet results in a rigid-body response for all the signals 
chosen. This maneuver was analyzed to obtain control derivatives for the rigid- 
body response described by Bqs. (15) to (31). The resulting match between the 
measured and computed response is shown in Fig. 42. The estimated control deri- 
vatives are in good agreement with those obtained from the maneuvers. The unmo- 
deled structural dynamic modes are evident, but it is apparent that the modes 
will have little effect on the rigid-body control derivatives. The differences 
between the measured and computed rigid-body responses (the residuals) for the 
time close to when the jets were fired are shown in Fig. 43. The data shown 
here are for a sample interval of 0.006 sec. Some persistent structural 
ringing is shown for the two rates and the lateral acceleration. However, when 
a jet is fired, the increased structural response is evident. The structural 
coefficients can be extracted directly from the residual as they were for the 
example in the previous section. It appears that there may be some contamina- 
tion caused by the rigid-body response at the instant the jets fire. If so, 
this contamination can be eliminated in one of two ways: either analyze the por- 
tion of the maneuver a tenth of a second after the jet fires, or adapt the equa- 
tions of motion to include the structural dynamics in addition to the rigid 
dynamics. The structural dynamics depicted in Fig. 43 have not been analyzed, 
but the procedure is straightforward. The procedure used on this case was the 
same as that used on the example in the proceeding section. It is apparent, 
however, that more than one structural mode would need to be included in the 
model. 

All the analysis techniques discussed in this paper apply to the analysis of 
this space shuttle example. If state noise is included in the mathematical 
model, then the linear form of Bqs. (5) to (7) would be required. In general. 
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if the structural partial differential equation can be expressed in tne linear 
foras of Eqs. (5) to (7) (with or without state noise), the structural nodes can 
be analyzed readily with the MMLE 3 program (Ref, 17) in the tine domain. If 

the analyst prefers, the problem can be expressed in the linear constant coef- 
ficient form and analyzed in the frequency domain, as described in Ref. 12. The 
relative advantages and disadvantages of time-domain analysis as compared with 
frequency-domain analysis are also discussed in that reference. If the equa- 
tions are nonlinear, but in the form of Eqs. (1) to (3), then maximum likelihood 
estimates can be obtained in the time domain. 


CONCLUDING REMARKS 


The computed simple aircraft example showed the basics of minimization and 
tha general concepts of cost functions themselves. In addition, the example 
demonstrated the advantage of low measurement noise, multiple estimates at a 
given condition, and the Cramer-Rao bounds, and the quality of the match between 
the measured and computed data. The flight data showed that many of these con- 
cepts still hold true even though the dimensionality of the cost function makes 
it impossible to plot or visualize. In addition, the techniques used for the 
aircraft problem were shown to be applicable to the flexible structure problem. 
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Table 1 Values of confuted time 
history with no measurement noise 


1 

6, deg 

P# deg/sec 

1 

0 

0 

2 

1 

0.9754115099857 

3 

1 

2.878663149266 

4 

1 

4.689092110779 

5 

1 

6.411225409939 

6 

1 

8.049369277012 

7 

1 

9.607619924937 

8 

0 

10.11446228200 

9 

0 

9.621174135646 

10 

0 

9.151943936071 


Table 2 Pertinent values as a function of iteration 


L 

Lp(D 

Lfi(L) 

♦ (L) 

♦ (L) 

J L 


0 

-0.5000 

15.00 

0.9048 

2.855 

21.21 


1 

-0.30C5 

9.888 

0.9417 

1.919 

0.5191 


2 

-0.2475 

9.996 

0.9517 

1.951 

5.083 x 

10” 4 

3 

-0.2500 

10.00 

0.9512 

1.951 

1.540 x 

10”® 

4 

-0.2500 

10.00 

0.9512 

1.951 

1.060 x 

10” 14 


Table 3 Values of computed time his- 
tory with added measurement noise 


1 

6, deg 

p, deg/sec 

1 

0 

0 

2 

1 

C. 4875521781881 

3 

1 

3.238763570696 

4 

1 

3.429117357944 

5 

1 

6.286297353361 

6 

1 

6.953798550097 

7 

1 

10.80572930119 

8 

0 

9.739367269447 

9 

0 

9.788844525490 

*0 

0 

7.382568353168 
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Table 4 Pertinent values as a function of iteration 


L 

Lp(L) 

bfiU) 

♦a> 

+(L) 


0 

-0.5000 

15.00 

0.9048 

2.855 

30.22 

1 

-0.3842 

10.16 

0.9260 

1.956 

3.497 

2 

-0.3518 

10.23 

0.9321 

1.976 

3.316 

3 

-0.3543 

10.25 

0.9316 

1.978 

3.316 

4 

-0.3542 

10.24 

0.9316 

1.978 

3.316 

5 

-0.3542 

10.24 

0.9316 

1.978 

3.316 


Table 5 

Mean and standard deviations for 

estimate* o* Lp 

Number 

of 

Sample mean. 

Sample standard 

Sample standard 
derivation of the 

cases , 

N 

U(Lp) 

deviation, o{Lp) 

mean, c(Lp)//N~ 

5 


-0.2668 

0.0739 

0.0336 

10 


-0.2511 

0.0620 

0.0196 

20 


-0.2452 

C.0578 

0.0129 


nble 6 Estimate of Lp and Cremer-Rao bound as 
a function of the square root of noise pouar 


Square root of 
noise power 

Estimate 
of Lp 

Cramer-Rao 

bound 

0.0 

-0.2500 


0.01 

-0.2507 

0.00054 

0.05 

-0.2535 

0.00271 

0.10 

-0.2570 

0.00543 

0.2 

-0.2641 

0.0109 

0.4 

-0.2783 

0.0220 

0.8 

-0.3071 

0.0457 

1.0 

-0.3218 

0.0579 

2.0 

-0.3975 

0. 1248 

5.0 

-0.6519 

0.3980 

10.0 

-1.195 

1.279 
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Roll rat# (p) 

Rolling momant (L # ) 



( c ) Simplified aircraft nomenclature * 
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Fig* 2 Concluded • 


Fig* 3 Time his- 
tory with no meas- 
urement noise • 
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Fig* 4 Comparison of measured and 
computed data for each of the first 
three iterations . 
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Fig* 5 rime his- 
tory vith meas- 
urement noise . 
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P* 10 
dag/aac 
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0 1 2 0 1 2 
Tima, aac Tima, aac 


Fig. 6 Comparison of measured and 
computed data for each iteration • 


20 


Roil rata, 

P, 10 
dag/aac 


0 1 2 
Tima, aac 

Fig* 7 comparison 
of estimated roll 
rate from no-noise 
and noisy cases • 





Fig. £ Cost function as a 

function if L p for no- noise case . 



Damping In roll, 

Fig. 9 Cost function as a function 
of Lp for noisy case • 
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Noiay 



Fig. 10 Comparison of the cost 
functions for the no-noise and noisy 
cases. 



Fig. 12 Comparison of newton- 
Raphson and Causs-Sf values of 
the second gradient the noisy 
case. 


Fig. 11 Gradient of J(L p ) as a 
function of L p for noisy case. 



Fig. 13 Comparison of first iter- 
ation step size for the newton- 
Raphson and Gauss-newton algorithms 
for the noisy c a»e. 
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Cmtmi mkm 
m4 mtknkmm 


Fig . I 7 M#tr of cost £um* 

tiom surface for fi>*j mime case* 
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Fif. J9 Detailed view of ox 
tion surface fox misy case 


Fig. 18 isoclims of cotwtarA cost 
in L n aM La fox the misy case- 




Pig. 20 isoclines of constant cost 
for region near minimum for noisy 
case. 



Fig. 21 isoclines and uncertainty 
ellipsoid of the cost function for 
the noisy case. 
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Fig. 23 Lateral -directional apace 
shuttle maneuver at a Mach number of 21. 
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Fig. 24 Examples of obtaining Lyj by 
simple calculations for the shuttle 
data from Fig. 23. 
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(c; Two iterations. 


(d ) Three iterations, 


Fig . 27 Continued. 
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Pig . 28 Variations of C n ^ with angle 
of attack without uncertainty levels. 


(e) Four iterations. 
Fig. 27 Concluded. 
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rig. 29 Variations of C Dp with angle of attack 
with uncertainty levels. 
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Fig. 30 Match of flight data obtained 
in turbulence (state noise ) and com- 
puted data obtained from maximum 
likelihood estimator that does not 
account for turbulence. 


Fig. 31 Match of flight data obtained 
in turbulence (state noise ) and com- 
puted data obtained from a maximum 
likelihood estimator that accounts 
for turbulence. 
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A 

Fig. 32 Power spectral density of $g 
obtained from maneuver shown in Fig . 
31. 



Frequency, Hz 

Fig • 33 Power spectral density of pg 
obtained from maneuver shown in Fig. 
31. 


Flight 



Fig • 34 Structural mode oscillation 
observed on the lateral acceleration. 



Fig. 36 Cost functional for 
10 cycles of data as function 
of frequency , showing close 
proximity of local minima 
to global minimum. 
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Fig. 35 Hatch of measured and com- 
puted lateral acceleration obtained 
when maximum likelihood estimator 
converged to local minimum . 
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Fig. 37 Simple scalar example illus- 
trating a local minimum similar to 
that shown for flight data in Fig. 35. 
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Fig. JS Simple scalar example showing 
only the first cycle . 



Fig. ?9 Cost function for one cycle 
of data as function of f zequency , 
showing wide region of convergence 
for global minimum • 
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Pig. 40 Acceptable match of Pleasured 
a net computed lateral acceleration . 
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Fig. 41 Dynamic response of space 
shuttle to firing of roll and yaw 
jets at an altitude of 430,000 Ft. 



TU.ia, aar 

Fig. 42 Maximum likelihood match 
of rigid-body response of the space 
shuttle . 
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Fig . 4J Differs 7 " , > en measured and computed rigid~body response 
( residual ) for «. shuttle. Altitude * 430,000 ft; dynamic 

pressure * 0* 
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